home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / diff / diff.c next >
Encoding:
C/C++ Source or Header  |  2001-03-15  |  4.7 KB  |  181 lines
  1. /* diff/diff.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 David Morrison
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_errno.h>
  24. #include <gsl/gsl_diff.h>
  25.  
  26. int
  27. gsl_diff_backward (const gsl_function * f,
  28.            double x, double *result, double *abserr)
  29. {
  30.   /* Construct a divided difference table with a fairly large step
  31.      size to get a very rough estimate of f''.  Use this to estimate
  32.      the step size which will minimize the error in calculating f'. */
  33.  
  34.   int i, k;
  35.   double h = GSL_SQRT_DBL_EPSILON;
  36.   double a[3], d[3], a2;
  37.  
  38.   /* Algorithm based on description on pg. 204 of Conte and de Boor
  39.      (CdB) - coefficients of Newton form of polynomial of degree 2. */
  40.  
  41.   for (i = 0; i < 3; i++)
  42.     {
  43.       a[i] = x + (i - 2.0) * h;
  44.       d[i] = GSL_FN_EVAL (f, a[i]);
  45.     }
  46.  
  47.   for (k = 1; k < 4; k++)
  48.     {
  49.       for (i = 0; i < 3 - k; i++)
  50.     {
  51.       d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i]);
  52.     }
  53.     }
  54.  
  55.   /* Adapt procedure described on pg. 282 of CdB to find best value of
  56.      step size. */
  57.  
  58.   a2 = fabs (d[0] + d[1] + d[2]);
  59.  
  60.   if (a2 < 100.0 * GSL_SQRT_DBL_EPSILON)
  61.     {
  62.       a2 = 100.0 * GSL_SQRT_DBL_EPSILON;
  63.     }
  64.  
  65.   h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2));
  66.  
  67.   if (h > 100.0 * GSL_SQRT_DBL_EPSILON)
  68.     {
  69.       h = 100.0 * GSL_SQRT_DBL_EPSILON;
  70.     }
  71.  
  72.   *result = (GSL_FN_EVAL (f, x) - GSL_FN_EVAL (f, x - h)) / h;
  73.   *abserr = fabs (10.0 * a2 * h);
  74.  
  75.   return GSL_SUCCESS;
  76. }
  77.  
  78. int
  79. gsl_diff_forward (const gsl_function * f,
  80.           double x, double *result, double *abserr)
  81. {
  82.   /* Construct a divided difference table with a fairly large step
  83.      size to get a very rough estimate of f''.  Use this to estimate
  84.      the step size which will minimize the error in calculating f'. */
  85.  
  86.   int i, k;
  87.   double h = GSL_SQRT_DBL_EPSILON;
  88.   double a[3], d[3], a2;
  89.  
  90.   /* Algorithm based on description on pg. 204 of Conte and de Boor
  91.      (CdB) - coefficients of Newton form of polynomial of degree 2. */
  92.  
  93.   for (i = 0; i < 3; i++)
  94.     {
  95.       a[i] = x + i * h;
  96.       d[i] = GSL_FN_EVAL (f, a[i]);
  97.     }
  98.  
  99.   for (k = 1; k < 4; k++)
  100.     {
  101.       for (i = 0; i < 3 - k; i++)
  102.     {
  103.       d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i]);
  104.     }
  105.     }
  106.  
  107.   /* Adapt procedure described on pg. 282 of CdB to find best value of
  108.      step size. */
  109.  
  110.   a2 = fabs (d[0] + d[1] + d[2]);
  111.  
  112.   if (a2 < 100.0 * GSL_SQRT_DBL_EPSILON)
  113.     {
  114.       a2 = 100.0 * GSL_SQRT_DBL_EPSILON;
  115.     }
  116.  
  117.   h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2));
  118.  
  119.   if (h > 100.0 * GSL_SQRT_DBL_EPSILON)
  120.     {
  121.       h = 100.0 * GSL_SQRT_DBL_EPSILON;
  122.     }
  123.  
  124.   *result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x)) / h;
  125.   *abserr = fabs (10.0 * a2 * h);
  126.  
  127.   return GSL_SUCCESS;
  128. }
  129.  
  130. int
  131. gsl_diff_central (const gsl_function * f,
  132.           double x, double *result, double *abserr)
  133. {
  134.   /* Construct a divided difference table with a fairly large step
  135.      size to get a very rough estimate of f'''.  Use this to estimate
  136.      the step size which will minimize the error in calculating f'. */
  137.  
  138.   int i, k;
  139.   double h = GSL_SQRT_DBL_EPSILON;
  140.   double a[4], d[4], a3;
  141.  
  142.   /* Algorithm based on description on pg. 204 of Conte and de Boor
  143.      (CdB) - coefficients of Newton form of polynomial of degree 3. */
  144.  
  145.   for (i = 0; i < 4; i++)
  146.     {
  147.       a[i] = x + (i - 2.0) * h;
  148.       d[i] = GSL_FN_EVAL (f, a[i]);
  149.     }
  150.  
  151.   for (k = 1; k < 5; k++)
  152.     {
  153.       for (i = 0; i < 4 - k; i++)
  154.     {
  155.       d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i]);
  156.     }
  157.     }
  158.  
  159.   /* Adapt procedure described on pg. 282 of CdB to find best
  160.      value of step size. */
  161.  
  162.   a3 = fabs (d[0] + d[1] + d[2] + d[3]);
  163.  
  164.   if (a3 < 100.0 * GSL_SQRT_DBL_EPSILON)
  165.     {
  166.       a3 = 100.0 * GSL_SQRT_DBL_EPSILON;
  167.     }
  168.  
  169.   h = pow (GSL_SQRT_DBL_EPSILON / (2.0 * a3), 1.0 / 3.0);
  170.  
  171.   if (h > 100.0 * GSL_SQRT_DBL_EPSILON)
  172.     {
  173.       h = 100.0 * GSL_SQRT_DBL_EPSILON;
  174.     }
  175.  
  176.   *result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x - h)) / (2.0 * h);
  177.   *abserr = fabs (100.0 * a3 * h * h);
  178.  
  179.   return GSL_SUCCESS;
  180. }
  181.